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Abstract 

We present a Monte Carlo algorithm that allows the simultaneous determination of a few ex- 
tremal eigenpairs of a very large matrix without the need to compute the inner product of two 
vectors or store all the components of any one vector. The new algorithm, a Monte Carlo imple- 
mentation of a deterministic one we recently benchmarked, is an extension of the power method. 
In the implementation presented, we used a basic Monte Carlo splitting and termination method 
called the comb, incorporated the weight cancellation method of Arnow et al, and exploited a new 
sampling method, the sewing method, that does a large state space sampling as a succession of 
small state space samplings. We illustrate the effectiveness of the algorithm by its determination of 
the two largest eigenvalues of the transfer matrices for variously-sized two-dimensional, zero held 
Ising models. While very likely useful for other transfer matrix problems, the algorithm is however 
quite general and should hnd application to a larger variety of problems requiring a few dominant 
eigenvalues of a matrix. 
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I. INTRODUCTION 



A common problem in computational physics is computing the eigenpairs of large matri- 
ces. We will present a new Monte Carlo algorithm that allows the simultaneous determina- 
tion of a few extremal eigenpairs of a very large matrix without the need to orthogonalize 
pairs of vectors to each other or store all the components of any vector. This algorithm is a 
Monte Carlo implementation of deterministic one we recently benchmarked [l|. The newly 
benchmarked algorithm is based on a refinement of the power method recently developed 
by Booth and does not require, as does the standard Ritz estimator , the ex- 

plicit computation of the inner product of two vectors. This feature is of special importance 
for Monte Carlo use because Monte Carlo sampling provides only successive estimates of 
eigenvectors represented by a very small subsets of their possible components. For such a 
situation, the explicit computation of an inner product is impossible. 

The basic power method is the traditional starting point for a Monte Carlo determination 
of the eigenpair associated with the eigenvalues of largest absolute value Ai. While various 
versions of the Monte Carlo power method often compute this dominant eigenvalue very 
well, computing subdominant eigenvalues X 2 , . . . has often proven much more difficult 
and is much less frequently attempted. Our Monte Carlo power method computes multiple 
extremal eigenpairs simultaneously. The particular algorithmpresented uses a basic Monte 
Carlo splitting and termination technique called the comb [g, |7[], incorporates the weight 
cancellation method of Arnow et al. [8], and exploits a new sampling method, the sewing 
method that does a large state space sampling as a succession of small state space 
samplings. 

In refining the power method we were targeting its use on matrices so large that they 
are unassailable deterministically because no single vector can be stored in memory. As the 
system size increases, finding a few extremal eigenpairs of the transfer matrix of the two- 
dimensional Ising model becomes such a problem. We will illustrate the effectiveness of the 
algorithm by determining the two largest eigenvalues of the transfer matrices for variously- 



sized two-dimensional zero field Ising models, exploiting Onsager's exact results [10|, LL3J for 
their values as benchmarks. We comment that two extremal eigenvalues of this matrix are 
of significant physical interest: the logarithm of Ai is proportional to the free energy, and 
the logarithm of the ratio A2/A1 is proportional to the reciprocal of the correlation length 
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that controls long range spin correlations near the critical point. Although our algorithm is 
extendable to finding more than two extremal eigenpairs, we will focus on finding just Ai 
and A 2 [lO, U\. 



In the next section, Section II, we summarize the basic features of the transfer matrix 
of the two-dimensional Ising model in a zero magnetic field. In subsequent sections we will 
reference these features to make our algorithm presentation more concrete. In Section III, we 
summarize our extension of the power method for the determination of the two eigenpairs 
corresponding to the two eigenvalues of largest absolute value. Then, in Section IV, we 
discuss the basics of our Monte Carlo implementation of this algorithm. We first use the 
algorithm for medium-sized matrices for which computer memory is adequate to store the 
eigenvectors and then use it for much larger- sized matrices for which it is not. For the 
Monte Carlo sampling of states in latter case, we use the sewing algorithm [9] that facilitates 
sampling of states from a large space from a smaller space. Results presented are for the 
determination of the two largest eigenvalues of the Ising transfer matrix for various lattice 
sizes. In the last section, we summarize our work and comment on its likely application to 
other systems. 

We note that the intent of the present application is presenting and benchmarking a 
new and relatively general numerical method and not numerically studying the finite-size 
scaling of the eigenvalues of the matrix used for the benchmarking. Because we reproduce 
the exact eigenvalues to satisfactory accuracy, our estimates will enjoy the same scaling as 
the well known exact results for this problem. When the magnetic field is not zero and 
the eigenvalues are not exactly known, various researchers have calculated up to the four 
largest eigenvalues by a deterministic version of the power method and have made extensive 



studies of the scaling of these eigenvalues 



12| . With the algorithm to be presented we have 



reproduced the first two of these eigenvalues and hence would reproduce the basic scaling. A 
more extensive study of the eigenvalues of the field dependent Ising model will be presented 
elsewhere 13| . The Monte Carlo approach increases the lattice sizes accessible. We note 
that these deterministic calculations required multiple processors. All our calculations were 
done on a single processor. 
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II. ISING MODEL 



In his much celebrated work Onsager calculated many of the properties of the two- 
dimensional Ising model exactly [h]]. His calculations started with the expression of the 
partition function in terms of its transfer matrix 
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151 ] . He then found all the eigen- 



values of this matrix analytically and showed that the scaling the dominant eigenvalue with 
the area in the thermodynamic limit (the area of the model approaching infinity) implied 
the onset of long-range ordering among the spin variables. The ratio of the second largest 
eigenvalue to the first is associated with the spatial behavior of this ordering. At the critical 
temperature, this ratio approaches unity as the area approaches infinity. 

We will consider an m x n Ising model defined with periodic boundary conditions in one 
direction and open boundary conditions in the other. The two-dimensional Ising model's 
energy is 

m— 1 n m n 

e {/4 = -j Vij^i+ij - J J2J2 Vij^ij+i C 1 ) 

i=l j=l i=l j=l 

Here, (i, j) are the coordinates of a lattice site. The Ising spin variable fit j on each site has 
the value of ±1, the exchange constant J > 0, and Hi )m +i = The symbol 

a j = (^l,j,^2,j, • • • , Mmj) (2) 

denotes a column configuration of Ising spins and there are 2 m possible configurations for 
each column. 



The transfer matrix A(a, a') follows from a re-expression of the partition function llj 



Z(m,m) = ^exp[-z/£({^})] 
M 

= A i (J i^ (7 2)A(a2,a 3 ) ■ ■ ■A{a m - U a m )A{a m ,a 1 ) 

CTl,...,(7 m 

= E^k^) (3) 

where v = J/ksT, ks is Boltzmann's constant, and T is the temperature, and A(a, a 1 ) is a 
2 m x 2 m matrix whose elements are 

(m— 1 \ / m \ 

v J2 Wk+iJ exp \ v J2 ( 4 ) 

As customary for Ising model simulations, we represent a configuration a by the first m bits 
of an integer between and 2 m — 1, with a set bit being a +1 Ising spin and an unset bit 



being a —1 spin. This convention maps the matrix element A (a, a') between configurations 
to an element Ay between integers. 

We note that A(a, a') is asymmetric and its elements are greater than zero so the matrix 
is maximally dense and hence irreducible. Because of the positivity and irreducibility the 
Perron-Frobenius Theorem [3] says that dominant eigenvalue is real, positive, and non- 
degenerate and all components of the corresponding eigenstate are real and have the same 
sign. The two largest eigenvalues of A, for finite m, are [l^] 



Ai = (2 sinh 2v) m l 2 exp 
A2 = (2 sinh 2z/) m//2 exp 



1 



(r)i + m H h m 



m— 1 j 



{r}2 + V* H 1" V2r. 



where 



irk 

cosh r/k = cosh 2u coth lv — cos — 



The transfer matrix of the Ising model can be symmetrized 
advantage for using this form. 



(5) 



(6) 



111 ] but we saw no computational 



III. POWER METHOD 



For some real-valued M x M matrix A, not necessarily symmetric, we will be concerned 
with the M eigenpairs (Aq,, ip a ) satisfying 



A^) a = X a ijj 



(7) 



In the simplest application of the power method [4j, an iteration is started with some ip, nor- 
malized in a convenient, but otherwise relatively arbitrary, manner and consists of iterating 
the two steps 

= Aip 
V = 0/||0|| 

until some convergence criterion is met. If we write 

M 



(8) 



(9) 



a=l 



and if |Ai| > | A2 1 > [A3 j > • • • > |Ajv|, then after n iterations 

M 



A n i) = Aj 



r / A \ n 

a=2 \ A 1/ 



(10) 
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Accordingly, as n — > oo, 



if) -> ^i/H^ill 

- Ai (11) 



Thus, the dominant eigenpair is simultaneously determined. For the norm of the vector cf) 
whose components are <pi, a frequent choice is 

II || = || ||oo= max (12) 

For deterministic calculations of a few dominant eigenpairs, say N, one of two approaches 
are typically tried. One approach is to use the power method to determine the dominant 
eigenpair, use deflation to project out this state out of the matrix, and then reuse the power 
method on the deflated matrix. To determine several eigenpairs simultaneously, the power 
method can be generalized to 

$ = Aty (13) 

where $ and $ are M x N matrices whose columns are orthogonalized to each other. 
This orthogonality needs maintenance throughout the computation or else all N vectors, 
represented by the columns of the initial ^, will converge to the one associated with the 
dominant eigenvalue j^]. This approach is called orthogonal [5] or simultaneous iteration 



I J 

For Monte Carlo calculations of a few dominant eigenpairs, we are proposing a quite 
different approach based on Booth's proposed refinement of the power method 0, 3]. This 
refinement uses the observation that for any eigenpair (A, tp) and for each non-zero component 
of the eigenvector, the eigenvalue equation Aif) = \ip can be rewritten as 

E AijTpj 

A = (14) 

and that similar equations can also be written for any number of groupings of components, 

E EAj^Pj E E •ly-'j E EAiV'i 

E A E A E A 

where the Ri are rules for different groupings. The groupings are quite flexible: they can 
overlap and their union need not cover the entire space. In addition, any two groupings, say 
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1 and 2, imply 

E^EE =E^«EE ^ (is) 

From L groupings of the components, Booth constructed L estimators for the L largest eigen- 
values and forced them to become equal by adjusting certain parameters at each iteration 
step. For the dominant two eigenvalues, we will do something similar. 

First we note that for almost any starting point ip = J^a^aipa, the power method will 
converge to (Ai,^i). To find the two extremal eigenvalues, we need two normalized, but 
not necessarily orthogonal, starting points ip' = J2 a u' a ip a and ip" = Ea^o^a jlj. At each 
step, we will apply A to them individually. Without any intervention both will project the 
same dominant eigenfunction so at each step we adjust the relationship between their sum 
to direct one to the dominant state and the other to the next dominant one. 

Formally, we start the iteration with ip — ip' + rpp" . Suppose at the n th step, xjj' and 
ip" have iterated to ip' and ip", then at the {n + l) th step we invoke (TT5T) and and (|T6il and 
require that 



E E Aijip'j + rj E E A^ E E A^- + v E E Ajf- 

teRi j «6-Ri j _ idR.2 j 1&R2 j 



(17) 



which leads to 



with 



Q2V +qiV + Qo = (18) 



<& = E € E E A^i - E € E E AhW 

idR2 i&Ri j ieRi idR2 j 

91= E#EE ±$ - E ft E E ^ 
+ Ei'EEM-Ei'EEAi* 

ieR,2 i&Ri j ieRi ieR,2 j 

* = E^EE a^ - E # E E (is) 

The algorithm is to apply A repeatedly until two real solutions rji and 7] 2 for (|18|) exist. One 
solution will be then used to guide further iterations to (Xi,tpi); the other, to (A 2 ,"i/>2)- 

The power method becomes [l|: choose two starting points ip' and ip", which need not 
be orthogonal, then for each iteration step, compute 

r - (20) 



and if the roots of (fTEj) are real [19j], update using 



iff' «- A0" + 7ftA.V>' 

^" <- A^" + r/ 2 ^ / (21) 



otherwise use 



«- At/>" (22) 



Eigenvalues can be estimated from 



E E A^Vj + Vi E EA^' 

E EA^ + 772 E E Art" 
Aa " E^ + ^E^ (23) 

ieRi ieRi 

where rji and ^2 generate the largest and next largest eigenvalue estimates. 



IV. MONTE CARLO IMPLEMENTATION 



The basic operation of a power method is a matrix-vector multiplication. Here, we 
now describe how we used the Monte Carlo method to estimate the repetition of such 
multiplications. 

In the basis defining the matrix elements of A, we write 

V = E^IO 

i 

< = £^'N> (24) 

i 

We will call the amplitudes io\ and uj" weights even though they are not necessarily all 
positive nor are the sums of their absolute values unity. We will also assume that the 
elements A^ of the M x M matrix A are easily generated, as we are ultimately interested in 
cases where M is so large that this matrix must be generated on-the-fly as opposed to being 
stored. Next, we imagine we have iV particles distributed over the M basis states defining 
A. Generally, iV <C M. Then, at each iteration step, we interpret Ay as the weight of a 



8 



particle arriving in state \i) per unit weight of a particle in state The action of Aij on a 
ip thus causes all particles currently in state \j) to jump to carrying to \i) their current 
weight u)j, modified by A^. 

The jumps will be executed probabilistically. To do this we let the total weight leaving 
state \j) be 

Wi = J>« (25) 

i 

and define the transition probability from \j) to \i) be 

Tij = Ay/Wj (26) 

The number Wj is called the state weight multiplier. How we use these densities will depend 
on the size and types of the matrices under consideration. 



A. Medium Matrices 



If M is sufficiently small so we can store all components of our vectors, then a Monte 
Carlo procedure for jumping is easily constructed. Instead of always (i.e., with probability 
1) moving weight A^ from state \j) to state |z), we will instead sample a |z) from Tj and 
multiply the transferred weight by the ratio of the true probability (1.0) to the sampled 
probability (Tij); that is, if state \i) is sampled, the weight arriving in state \i) from \j) is 
multiplied by 

^r^^kr™' (27) 

As for many Monte Carlo simulations, as is the case for the transfer matrix of the Ising 
model, the particle weights defining the eigenvector associated with the largest eigenvalue 
can be made all positive. The second eigenfunction however must be represented by some 
particles of negative weight and some particles of positive weight. These negative and posi- 
tive weights must for some jumps at least partially cancel to maintain a correct estimation 
of the second eigenfunction. When N <C M, this cancellation does not occur often enough 
in a Monte Carlo simulation without proper design: because the number of states vastly out 
numbers the number of particles, the probability that a negatively and a positively weighted 
particle randomly arrive in the same state becomes trivially small. 

There are several ways to arrange the cancellation jjj. We found the Arnow et al. [8] 
algorithm effective and convenient. First, we consider two particles of weights Wi and W2 in 



states \ji) and | j'2) , and then let and Ty 2 be the probabilities of their reaching state |z). 
Next we use the weight multiplier Wj for jumping from state \j) to state \i) and suppose 
that states and ^2) are sampled from the density T i j 1 + Tij 2 . This density can be sampled 
by sampling a from T i j 1 and a ^2) from T^ 2 . The true probability that particle 1 jumps 
to \ii) is T il j 1 and the true probability that particle 2 jumps to is T il j 2 . The ratio of the 
true density to the sampled density for particle 1 is 



T ■ 



T 4- T ■ 

so that the weight arriving at state from particle 1 is 

wiW h T ilh 



(28) 



(29) 



The true probability that particle 2 arrives at |zi) is T il j 2 . The ratio of the true density to 
the sampled density for particle 2 is 



T ■ 

A tin 



so that the weight arriving at from particle 2 is 

w 2 W h T ilj2 



T 4- T ■ 

Thus the total weight arriving at is 

WiW h T ilh + w 2 W j2 T ilj2 
T 4- T ■ 

By similar arguments, the total weight arriving at ^2) is 

w x W h T i2h + w 2 W j2 T i2j2 



(30) 



(31) 



(32) 



(33) 



We note that to get meaningful cancellation, say between particle 1 with weight w± < 
and particle 2 with weight w 2 > 0, the transition probabilities must overlap somewhat. For 
example, if in (I3"3"|) Tj 2J1 = e<l, then the total weight arriving at i 2 is essentially just the 
weight arriving of particle 2 alone. 



w{W h t 4- w 2 Wj 2 T i2j2 



£ 4- T{ 2 j 2 



w 2 W j2 (34) 



How one arranges better overlap is problem dependent. For our Ising simulations, we sorted 
the particles into state order (a state is represented by the bits of an integer). Particles 1 
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and 2 were then sampled together according to the Arnow et al. scheme, then particles 3 
and 4 are sampled together, and so forth. The fact that the list is ordered means that there 
are (typically) many nearby states \i) that are accessible from both particles i and i + 1 
with nontrivial transition probabilities. 

As the iteration progresses, the absolute value of the weights of some particles becomes 
very large, and those of some others, very small. As standard for Monte Carlo methods with 
weighted particles, particles with weights of small magnitude are stochastically eliminated 
and those with large magnitudes are stochastically split. To do this we used a procedure 
called the comb [f],Q]. It is described in Appendix A. 

The steps of the algorithm are: First, we initialize the states and weights of two vectors. 
For the Ising simulation, each vector had the same states but different weights. We selected 
the states uniformly and randomly over the interval (0, 2 m — 1) and selected the io[ uniformly 
and randomly over the interval (0,1) and the to" uniformly and randomly over the interval 
(-0.5,0.5). Then, for a fixed number of times we iterate. At each iteration we execute the 
jump procedure for each particle, place the particle list in state order, effect cancellations, 
estimate the eigenvalues from ( 1231) . update if)' and ip", and then comb. R\ consisted of the 
states for which more than half of its m bits were 0's and R2 consisted of the states for 
which more than half of its m bits were l's. Additional algorithmic details are given in 
Appendix B. 

Table [J shows the computed Ai and A2 for m = 12. A computer program was written 
to make 20 independent (different random number seeds) calculations using N particles per 
iteration, for various values of N. Note that even when N = 100 << M = 2 12 = 4096 
the method still separates the eigenfunctions. The simulations were run at the critical 



temperature of the infinite lattice; that is, v = 0.4406867935097715 [TJJ. The last line in the 
table gives the eigenvalues deterministically obtained by using ((Sj). From the last few lines, 
we also see that when the number of particles floods the number of states, that is, when all 
the basis states are being used multiple times, exceptional accuracy is obtained. 

Also in the table is the timing for each run. The runs were done on a single 1.5 GHz Mac 
PowerPC processor. We note that the runs appear to scale sub-linearly with the number of 
particles. This scaling is deceiving. For a small number of particles, the time required to 
set up the transition matrix is significant compared to the rest of the calculation. The 
Monte Carlo part is dominated by the state order sort which scales as iVlogiV, so as N 
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TABLE I: Particle number N dependence of individual runs for the m = 12 Ising Model. The Aj 
and Gi {% = 1,2) are run's eigenvalue averages and error estimates. The row labeled "Onsager" are 
the eigenvalues obtained from ([5]). Also given is the wall clock time for each run. The order of the 
transfer matrix is 4096. 



N 


Ai 


0"l 


A 2 


°2 


minutes 


100 


71415.164 


65 


67021.909 


103 


5.4 


1000 


71527.110 


17 


67023.420 


31 


7.9 


10000 


71553.325 


5.3 


66956.314 


9.1 


22.8 


100000 


71557.854 


2.0 


67005.486 


3.2 


150.0 


1000000 


71557.129 


0.36 


67010.989 


0.58 


1474.0 


Onsager 


71557.047 




67010.869 







becomes very large, the run time scaling should eventually be slightly super-linear. 



B. Large Matrices 

For M < 2 12 , sampling from the cumulative probability C, = X)fc=o^*j works well. If the 
number of states gets too large (M = 2 12 was our limit), then C, cannot be sampled directly 
because it cannot fit in the computer's memory. In this case, we could just randomly pick 
from any state \j) any state \i) with probability 1/M instead of always picking a state \i) 
(i.e., picking it with probability 1). Then, if state \i) is sampled, the weight of particles 
arriving in state \i) is 

Aa—^r = AaM (35) 
3 1/M 3 v ; 

The problem with this approach is that the Ay can have immense variation so that this 
simple sampling scheme is unlikely to work well as a Monte Carlo method. This situation 
is especially true for the Ising problem. A large part of such variations however can be 
removed by sampling the new state in stages and then sewing the stages together. 

To explain the sewing procedure [9(, we will first assume that we can write any state \i) 
in our basis as a direct product of the states in a smaller basis, for example, 

\i) = \i2) In) 



12 



Instead of transferring weight 



W* = 5>y (36) 

k 



from state \j) to state \i) with probability 

T V = A ij / y £ l A kj = A 1j /W j , (37) 
fc 

we will use the that would apply to the smaller set of states and then make an appropriate 
weight correction. 

For each smaller set of states, we rewrite the analogous transition probability from state 
\j) to state \i) as 

t ij = a ij /w j , (38) 

and the analogous weight multiplier as 

w i = £ a kj (39) 
k 

We thus will sample and ^2) from the probability function 

Now, we define CV,- to be the weight correction necessary to preserve the expected weight 
transfer from state \j) = \j2)\ji) to state It satisfies 

Aij = Cijti 2 j 2 ti 1 j 1 (41) 

Using ([38]) in dUD, we find 

A = c a ^h^n (42) 

Thus 

A ■ 

C ij = w jl w h - l j— (43) 



This sewing method generalizes easily. For k sets of states, (jlzj) and (|43j) become 

4* = Cy f[ (44) 
n=l 

and 

fe it) • 

a, = a^ n (45) 



Tl=l ^^njn 
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For the Ising problem, the weight correction is 

k 

dj = exp(z/A) II w in (46) 

n=l 

where vDi is the energy difference per fcgT between calculating with the bits together and 
the bits separately [9|. 

TABLE II: Estimates of the two dominant eigenvalues Ai and A2 and their errors of the transfer 
matrix for variously-sized two-dimensional Ising models. Each estimate was based on 20 indepen- 



dent simu 
result [1CJ, 



ations. Also given for each lattice size is value for Ai computed via Onsager's exact 

li. 



Matrix Size 


Ai (Onsaj 


;er) 


Ai 


Ai (Onsager) 


A 2 


2 i6 x 2 16 


2.93297 x 


10 6 


2.93307 ± 0.00008 x 10 6 


2.79225 x 10 6 


2.79482 ± 0.00010 x 10 6 


2 24 x 2 24 


4.95473 x 


10 9 


4.95480 ± 0.00020 x 10 9 


4.79510 x 10 9 


4.79502 ± 0.00029 x 10 9 


2 32 x 2 32 


8.39316 x 


10 12 


8.39311 ± 0.00049 x 10 12 


8.18959 x 10 13 


8.18807 ± 0.00061 x 10 12 


2 40 x 2 40 


1.42333 x 


10 16 


1.42334 ± 0.00007 x 10 16 


1.39565 x 10 17 


1.39558 ± 0.00008 x 10 16 


2 48 x 2 48 


2.41504 x 


10 19 


2.41522 ± 0.00019 x 10 19 


2.37584 x 10 19 


2.37481 ± 0.00054 x 10 19 



Using this sewing algorithm for the sampling of states, we computed the first and second 
eigenvalues and their standard deviations for variously sized two-dimensional Ising models 
by sewing sets of 8 bits. The results are shown in Table [III Twenty independent runs 
were done for each size with 500 iterations per run. There were 1 million particles per 
iteration for m = 16, 24, and 32 and 5 million for m = 40 and 48. Only the second half of 
the iterations for each run used in the estimation process. This choice for "burn-in" was 
arbitrary and excessive. Typically the iteration converges to its fixed point in about 10 or 
fewer steps. Clearly, the estimated mean of the eigenvalues are consistent with the analytic 
result of Onsager. The m = 48 run took about 2.15 hours on a single 1.25 GHz Alpha EV6 
processor, v = 0.4406867935097715, the bulk critical value. R± consisted of the states for 
which more than half of its m bits were 0's and R2 consisted of the states for which more 
than half of its m bits were l's. In Table ILTTj we present the run averages for the M = 2 48 
case. We note that 2 48 « 2.8 x 10 14 . 
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TABLE III: Run dependence of Ai and A 2 for the 48 x 48 Ising Model. Each run had 5 million 
particles and 500 iterations. Only the last 250 were used to computed the stated run averages. 



Run 


Ai 






A 2 






1 


2.415237261676499 


X 


10 19 


2.376562826996199 


X 


10 19 


2 


2.416053055304024 


X 


10 19 


2.376357592801913 


X 


10 19 


3 


2.414208101931427 


X 


10 19 


2.373861137045189 


X 


10 19 


4 


2.414240257597889 


X 


1Q 19 


2.374603969362707 


X 


1Q 19 


5 


2.415129736468956 


X 


10 19 


2.375668969910789 


X 


1Q 19 


6 


2.414540994327886 


X 


1Q 19 


2.367761547315434 


X 


1Q 19 


7 


2.415931928295190 


X 


1Q 19 


2.371678054324882 


X 


1Q 19 


8 


2.416760807133448 


X 


1Q 19 


2.380200781153834 


X 


1Q 19 


9 


2.413931047178054 


X 


1Q 19 


2.374313310822579 


X 


1Q 19 


10 


2.414288017604171 


X 


1Q 19 


2.374183202851979 


X 


1Q 19 


11 


2.415907427261585 


X 


1Q 19 


2.374280853036757 


X 


1Q 19 


12 


2.415718175365685 


X 


1Q 19 


2.374828531467586 


X 


1Q 19 


13 


2.415613030894034 


X 


1Q 19 


2.374305559747647 


X 


1Q 19 


14 


2.414272502002316 


X 


1Q 19 


2.374687824232260 


X 


1Q 19 


15 


2.416636076568460 


X 


1Q 19 


2.376864186656784 


X 


1Q 19 


16 


2.414631286368427 


X 


1Q 19 


2.376097488856122 


X 


10 19 


17 


2.415202515765169 


X 


1Q 19 


2.376875902411582 


X 


1Q 19 


18 


2.415187019163818 


X 


10 19 


2.373320142080186 


X 


1Q 19 


19 


2.414692528984318 


X 


1Q 19 


2.374355064690523 


X 


1Q 19 


20 


2.416154446732387 


X 


1Q 19 


2.375485037469832 


X 


1Q 19 



V. CONCLUDING REMARKS 

We have presented a Monte Carlo algorithm that enables the determination of two ex- 
tremal eigenvalues of a very large matrix. The explicit demonstration of the power of our 
algorithm was the determination of the two largest eigenvalues of the transfer matrix of the 
two-dimensional Ising model in a zero magnetic field. The convenience of this matrix was 
the existence of exact expressions for its eigenvalues for any finite-sized system. We were 
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able to reproduce the exact values for the two dominant eigenvalues to within the statistical 
error of our simulations. 

Physics and chemistry provide numerous problems where obtaining several extremal 
eigenvalues is important. One problem class is finding the ground state and a few excited 
states of atoms, molecules, solids, and nuclei [20[. For these quantum problems, quantum 



21 



22] are commonly used. In fact, 



8] method used here was first proposed as 
21] . We comment that our new algorithm 



Monte Carlo Monte Carlo (QMC) projector methods 
it was for these methods that the Arnow et al. 
a method for taming the fermion sign problem 
requires sampling from the second eigenfuction which must have positive and negative com- 
ponents just as a fermion state must also have. Our successful cancellation of signed particles 
representing this eigenfunction and the slightly superlinear scaling of our computation time 
with system complexity, as opposed to an exponential scaling, demostrates that we have 



"solved" the sign problem |23j for our app 
We also note that several investigators 



ication. 



241 ] have adapted the technology of the diffusion 



QMC method 21] to find the dominant eigenvalue of the transfer matrix of various classical 
spin models. This Monte Carlo approach is distinctly different from ours even if ours were 
restricted to just the dominant state. For example, we do not use an importance function 
to guide the sampling nor do we use back propagation to enhance estimators. The current 



extension of diffusion QMC to the concurrent calculation of multiple eigenpairs [21|, [25| is 
also quite different from our extension. Diffusion Monte Carlo analogs the simultaneous 
iteration method mentioned in Section III. (Still another approach to concurrent eigenvalue 
estimation was used by Hasenbusch et al. 26].) We further note Nightingale and Blote's 



271 ] use of a Monte Carlo power method to find the second eigenvalue of a Markov chain 



transition matrix satisfying detailed balance. For this case the dominant eigenpair is known 
a -priori and was used to deflate the matrix so the projection was to the second largest 
eigenvalue. We comment that Booth has presented an algorithm for determining a single 
eigenvalue without the need to know or concurrently determine any other one. 

Our transfer matrix was positive, asymmetric, and dense. How is our algorithm changed if 
the matrix lacks one or more of these properties? For simple test cases, we have successfully 
constructed deterministic procedures for matrices whose elements are complex valued. Also 
for simple test cases, we have had success for real asymmetric matrices whose eigenvalues 
are complex valued. Devising Monte Carlo algorithms for real, symmetric, sparse matrices 
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has however received more of our attention 



28]. 



For an indefinite but symmetric matrix, the symmetry and reality assure the reality of 
the eigenpairs, but do not fix the sign of the dominant eigenvalue or the components of 
its eigenvector. Such matrices can also have the dominant eigenvalue being degenerate. In 
our algorithm a signature of the dominant eigenvalue being degenerate is the quantity q\ 
( TT9|) in our quadratic equation (fl8|) becoming extremely small. In cases where we knew the 
degeneracy a priori, we saw our estimate of Ai being very close to A2, that is, we saw our 
estimates of the dominant degenerate eigenvalue being slightly but artificially separated. 
This type of difficulty is inherent to a power method. If the matrix elements are of mixed 
sign, the modifications to the algorithm are quite simple: we take the absolute value of the 
and multiply the weight multiplier by the sign of A^. We also have to account for the 
sign of the weights uj[ and u'( appropriately. Weight cancellation must be done with care. If 
the matrix is sparse, other options for selecting the state to jump to become available that 
avoid storing the cumulative distribution, improve efficiency, and reduce variance 28]. 

Here, we have focused principally on the determination of eigenvalues, noting the rapid 
convergence to their values, but what about the determination of the eigenvectors? If we can 
store the vectors in computer memory we can accurately determine them. For a problem in 
a continuum, for example, finding the dominant eigenvalues of an integral equation such as 
those occurring for the transport of neutrons, we have accurately estimated the two dominant 
eigenfunctions with an efficiency enhanced over the basic power method for determining just 
the dominant eigenf unction 291 ] . 

The dynamics of determining eigenvectors is as follows: Because the standard power 
method for the dominant eigenpair converges to ip\ as A2/A1 ( TTOjl . the dominant eigenvector 
becomes increasingly difficult to determine as Ai and A2 become very close in magnitude. 
Because we are projecting to the first two dominant eigenpairs, our method converges to ipi 
as A3/A1. This is usually a gain relative to the standard power method. Convergence to ip2 
goes as A3/A2. This ratio illustrates the fact that in principle we can accelerate convergence 
to ipi by seeking the L > i highest eigenfunctions and thus converging to ipi as 
To do this, we need to use L groupings and L starting states and adjust the sum of the 
iterated states to guide them to convergence to L different eigenvalues. We have observed 
accelerated convergence for the L = 3 case on simple test cases. We comment that finding 
the eigenvectors of the transfer matrix of the Ising and other models is a problem for which 
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finding the eigenvectors becomes increasingly difficult as the system size increases, if we could 
in fact store them. For this type of matrix, the eigenvalue spectra is expected to become 
quasi-continuous starting with the second eigenvalue and the first and second eigenvalues 
become asymptotically degenerate at a critical point [30t ] . 

In closing, we believe that the algorithm presented here is accurate, easy to implement, 
and applicable to many other problems. Different problems afford the opportunity to im- 
prove its efficiency by modifying some of its details. Wider use of the algorithm will define 
more crisply its strengths and limitations than is possible by just the present application. 



APPENDIX A 



The comb is a stochastic procedure for selecting M particles, not necessarily all different, 
from a list of N unevenly weighted particles and preserving the total weight. If uoi > is 
the weight of an individual particle in the original list, its weight in the new list becomes 
Wt/M where Wt = Ef ^% is the total weight. To effect this procedure, we construct the 
cumulative sums of the original weights 

c j = j2 uj i ( A1 ) 

i=l 

where j = 0, 1, 2, . . . , N and Cq = and Cn = Wt, draw a random number £ uniformly 
distributed over (0,1), and then for k = 1, 2, ... , M, select particle j for the new list if 

Cj _ 1< {k-l+i)^<C j (A2) 

Dependent of the size of the difference Cj-\ — Cj relative to Wt/M, particle j is selected 
zero, one, or more times. The procedure calls the random number generator only once, and 
the new particle list is always of a predetermined size. 

In the present simulations we have a list of particles and to each particle % is associated a 
state \i) and two weights, uj[ and u" . If we were to comb the two lists of weights separately, 
we would in general produce two lists of states. We wanted one such list. Accordingly, we 
combed in the following manner: we formed 

Pi = Tr*L- (A3) 

E H\ 



18 



p'l = -P^- (A4) 

£ Kl 

i=l 

and then generated the cumulative distribution. 

C^tliPi+P") (A5) 

i=l 

Now we combed the Cj as before, noting that we have chosen a normalization such that 
Wt = 1. Instead of giving the selected particle j a weight Wy/M for both to'- and w", we 
instead assign 

^sign(^) 

w i = — ; ir~ 

3 ti+P" 

" (A7) 

APPENDIX B 

The following are some useful particulars of our Monte Carlo implementation of our 
modified power method. First, instead of updating via (1211) we update via 

i[>' <- Aip' + ^Aip" 

tfj" <- —Aifj' + Aifj" (Bl) 

V2 

This form is more symmetric and makes more explicit that we are trying to remove ip" from 
ip' and vice versa. To avoid over correcting we have sometimes found it useful to update via 

if)' <- Aip' + T&Aip" 

<- — A^' + Aip" (B2) 
m 

where 

7i[ = sign(?7 1 )min(a, |r]i|) (B3) 
r]' 2 = sign(?72) min(a, \r] 2 \) (B4) 

and a is some small positive number. 
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Besides the eigenvalue estimators (|23|) . there are others (fT5l) equally valid that provide 
multiple cross-checks to fl23l) . These include 



E E Aj^j E E Aji>j 

Ai= e« = ^*r (B5) 

E E.l,'.;J E E-l,<.;" 
A2= E< = Etf (B6) 

all of which involve sums needed for (123]) so they cost nothing extra to compute. Another 
use of these estimators is monitoring the utility of the regions R\ and R2, as an estimator 
will perform inconsistently if too few particles occupy a region or if the sum of the weights 
in a given region is small. The most useful regions are those that are a major subset of the 
positive and negative regions of the second eigenvector. 

For the Ising model we obtained some additional useful weight cancellation by exploiting 
the fact that most of the particles are in fully magnetized states and before using the comb 
replacing all the particles in the "up" state by a single particle whose weight is the sum of 
the "up" weights of these particles and doing similarly for the particles in the "down" state. 

Because the should be positive, after combing we explicitly set each uj[ equal to its 
absolute value. THe step is useful because even though we initialize ip' to have all positive 
components, the weight cancelation procedure can make some of them negative. Once so, 
they might project out rapidly and may adversely affect the eigenpair estimation until they 
do. 
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